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Abstract It has now become customary in the field of numerical relativity to couple 
high order finite difference schemes to mesh refinement algorithms. To this end, differ¬ 
ent modifications to the standard Bcrger-Oliger adaptive mesh refinement algorithm 
have been proposed. In this work we present a fourth order stable mesh refinement 
scheme with sub-cycling in time for numerical relativity. We do not use buffer zones to 
deal with refinement boundaries but explicitly specify boundary data for refined grids. 
We argue that the incompatibility of the standard mesh refinement algorithm with 
higher order Runge Kutta methods is a manifestation of order reduction phenomena, 
caused by inconsistent application of boundary data in the refined grids. Our scheme 
also addresses the problem of spurious reflections that are generated when propa¬ 
gating waves cross mesh refinement boundaries. We introduce a transition zone on 
refined levels within which the phase velocity of propagating modes is allowed to de¬ 
celerate in order to smoothly match the phase velocity of coarser grids. We apply the 
method to test problems involving propagating waves and show a significant reduction 
in spurious reflections. 

Keywords Numerical Relativity- Higher Order Mesh Refinement- Interface Boundary 
Conditions- Fixed Mesh Refinement 


1 Introduction 


Long term stable evolution of non linear hyperbolic partial differential equations of¬ 
ten require techniques to efficiently deal with vast length scales while simultaneously 
resolving fine scale features. Indeed, many numerical simulations in computational 
astrophysics, cosmology, numerical relativity and fluid dynamics are confronted with 
processes that span a wide range of time and length scales. In the context of numer¬ 
ical relativity, these simulations are often performed in full three spatial dimensions 
without symmetry assumptions. For these cases, running fine unigrid integrations is 
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computationally expensive and often impractical. Recent establishments in the numer¬ 
ical solution of partial differential equations by hnite difference techniques has seen 
an increasing use of nested grids and mesh refinement techniques in order to minimize 
the truncation error incurred with minimal computational and memory requirements 


The principle behind mesh refinement schemes is to recursively refine areas of the 
computational grid that are likely to induce higher discretization errors. This ap¬ 
proach efficiently focuses computational effort and resources in places where it is 
needed compared to refining the entire grid. Extensive theory has since been devel¬ 
oped for the method in different contexts, and mesh refinement algorithms have been 
widely adopted in the literature pnil23ll291ll2llM] . One of the key aspects in the im¬ 
plementation is the way in which the decision to add or remove levels of refinement is 
made. In this work we focus on the concept of fixed mesh refinement (FMR), where the 
grid hierarchy is created once and remains fixed for the duration of the computation 
mmm- This differs from adaptive mesh refinement (AMR), where the algorithm is 
endowed with error estimation routines that dynamically determines which areas need 
refinement. For most applications in numerical relativity, one can know before hand 
which areas within the computational domain require more refinement. As a result the 
FMR conhguration is widely adopted in the held. In the context of numerical relativ¬ 
ity, mesh refinement has seen application ranging from studies in critical phenomena, 
single black hole space times, binary black hole collisions, black hole lattice universes, 
neutron stars and core-collapse supernovae, within the non linear framework of general 
relativty, see for example pTl 1291 17] and references therein. 

Another component of the scheme is the inter-level coupling among nested grids. 
While the coarse grid has to supply boundary conditions to finer grids during evolu¬ 
tion, one can choose not to update the coarse grid solution with the fine grid solution. 
In this case, the coarse grid solution is independent of that of the nested grid, a con¬ 
hguration that is referecl to as one way (parasitic) nesting. In this work we employ 
two way (interactive) nesting, where we, in addition, update coarser levels with finer 
levels once the finer levels have been integrated to the same time level as the coarser 
ones. See [T7] for a comparison between parasitic and interactive coupling within the 
context of the shallow water equations. 

Traditionally mesh refinement techniques were coupled to second order convergent 
methods. On the other hand, recent trends in the numerical simulation community 
has seen the coupling of higher order hnite difference methods to the mesh refinement 
framework |30llT9lfTT[l36L[25] . This combines the efficiency of local mesh refinement with 
the robustness and accuracy of higher order methods. However, there is an inherent 
incompatibility between high order time discretization schemes with the standard 
mesh refinement algorithm that may result in loss of convergence or even instabilities 
[El- This issue is related to how the computation of boundary data for the refined 
grids is handled. A search for a stable high order mesh refinement implementation 
has resulted in several modifications to the standard method in an effort to address 
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this subject. For example, PSIE] employs the idea of buffer zones, where additional 
points are added closer to the boundary of the refined gridi-0 In this setup, boundary 
conditions are not prescribed explicitly in the refined levels, the integration is only 
applied to a progressively smaller domain in the refined grids and the buffer zone is 
ultimately discarded. Another approach is the tapered boundary approach mm- 
Here, one performs the integration at level l using the past domain of dependence of 
the child grid only. Other approaches have been to refine only in space and using the 
same time step for all levels [TO]. In this work we use a framework where we refine 
both in space and time and the treatment of interface boundaries is dictated by the 
time marching algorithm, fourth order accurate Runge Kutta algorithm in this case. 

In addition to issues of convergence and stability, one has to address the problem of 
spurious reflections off refinement boundaries that arise when propagating waves cross 
refinement boundaries. This is essential for gravitational wave source simulations as 
the waves are normally extracted at a large radius. As a result, propagating waves will 
have crossed several refinement boundaries, before reaching the radius of extraction. 
In [ID] , the idea of derivative matching was proposed in order to minimize spurious 
reflections for second order convergent schemes. In addition, the concept of mesh 
adapted stencils (MAD) was introduced in [2]. However these implementations do not 
involve refinement in time. Other methods that have been applied in the Advanced 
Weather Research and Forecasting Model, Euler equations and Maxwell equations 
involve the use of sponge layers in the refined levels to ensure that the solution in the 
refined levels will be nudged towards that of the coarser grids at refinement boundaries. 
This may involve the addition of artificial damping and dissipation terms in the system 
under consideration mm- See also the treatment of m in the case of first order 
convergent schemes. In this work we propose a simple scheme that is adopted from 
the animation and image processing community |14j . to deal with transitions from 
fine to coarse grid solutions at refinement boundaries. 

Finally, we note that high order mesh refinement implementations are often en¬ 
dowed with several lower order approximations in order to reduce the overall compu¬ 
tational demand [71 1251136] . In the case where there is refinement in time, this often 
includes the use of lower order interpolations in time. This results in some savings 
in memory usage as the storage of past time levels is lower for lower order interpo¬ 
lations. For example, one typically needs n + 1 past points in order to have a n-th 
order polynomial interpolant. In addition, one can reduce the number of buffer zones 
in the refinement boundaries by successively lowering the order of finite differencing 
as one approaches the refinement boundary. Other approaches in the case where ar¬ 
tificial dissipation is needed, include the use of Kreiss-Oliger dissipation operators of 
order lower than that required by the finite differencing operator. This lowers the 
number of required ghost points as well as being cheaper computationally. Ultimately, 
these approximations have a bearing on the convergence order, the amount of spuri- 

1 It is to be understood that the buffer zone is distinct from the ghost zone, which is traditionally used to hold boundary data 
in finite difference schemes. 
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ous reflections induced at interface boundaries and the overall performance of a given 
implementation. In order to keep the overall scheme consistent, we make no such 
simplifying assumptions in this work. 

This paper is organized as follows: in Section [2] we convey the framework of our 
FMR approach. We review our boundary application method in £ |2.4.1 and introduce 
the transition zone in 1 2.4.2 We present our results in Section [4] and finally concluding 
in Section 0 


2 Generalities 

2.1 Grid layout 


The grid hierarchy is arranged by first discretizing the spatial domain into a relatively 
coarse uniform mesh that covers the entire computational domain. This constitutes 
the base or root grid H° 0 with mesh size Jiq. Finer grid patches of mesh sizes hi/r are 
then overlaid as required to the base grid with each grid at level l having mesh size 
hi_i/r. More than one grid patch can be added in a given level. This forms a tree or a 
hierarchy of grids H l , where the indices l and p represent the level and patch number 
respectively. This configuration is depicted in Figure [T[ It is at this point that we 
emphasize that each grid in a given level l has its own solution vector and is evolved 
independently of all other grids. Of course it has to depend on the parent grid, within 
which it is nested, for boundary data. 

Each grid patch added at a given level must satisfy certain conditions. Among 
them, the idea of proper nesting: A fine grid at level l must start and end at the 
corner of a cell belonging to level l — 1. Moreover, grids at higher levels cannot ‘float’. 
This means that if there is a grid at level l + 2, it must be contained in a grid at level 
l +1 that is itself properly nested on a grid at level l. The refinement factor r must be 
an integer, and is the same for all levels. This results in a constant CFL for all added 
levels, thus the same integration routine of the base level is stable on all levels, if it 
is stable on level 0. It also implies that grids at higher levels require r l time steps to 
catch up with a single time step of the base grid. 


2.2 Inter-level communication 


Each grid in a given level l can be indexed independently by its own (p, T, ki ) coordi¬ 
nate system. However, for reasons of inter-level communications, there is a mapping 
from the H l p coordinates ( H,jhh ) to the H l ~ 1 coordinates (p_i,j;_i, h-i) and vice- 
versa. This relation is expressed as, 

H — mod(p, r ) m 

H- 1 = - 1-U 

r 

for a staggered hierarchy. Such communications are necessary for the computation of 
initial conditions, boundary conditions and for updating the coarse grid with the fine 
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Fig. 1 A grid hierarchy demonstrating proper nesting. A single grid H° 0 covering the entire domain is marked by ‘a ’. There 
is one refinement grid H 1 0 at level one marked by ’b’. Two disjoint refinement grids are H 2 0 and H 2 1 are marked by ‘c’ and 
‘d’ respectively. Note that the ghost zones are not included in the grids. This figure is used to emphasize that the grid hierarchy 
is not some complex data structure, but that the overlaid grids are independently stored in memory. 


grid solution. Initial data can be generated by spatial interpolation from the previous 
grid level ( prolongation ) or by calling the same initialization routine that was used to 
initialize the base grid. All levels are added and initialized at the same initial time 
to- Once all levels have been integrated to the same time, data in the finer meshes is 
used to update data in the coarse levels through the use of interpolation operators, a 
process called restriction. 

It is important to note that with the mapping of indices |TJ all fine grid points 
are staggered about coarse grid points at lower levels. This has the advantage that 
if the base grid is discretized strategically to ‘avoid’ certain points (. x,y,z ) € R 3 , 
(e.g. by using a cell centered grid to avoid dealing explicitly with the points on the 
edges) such points remain avoided in all refined levels. However, from a computational 
standpoint, this may be expensive because communicating data across levels requires 
three dimensional interpolation all the time. This is different from cases where some 
points in the fine grids are allowed to coincide with the ones from coarser levels. For 
such cases, inter-level communication can occur via injection , where data is simply 
copied to corresponding points at a given level. 


2.3 Numerical Integration 


Spatial differentiation is handled through fourth order convergent finite differencing. 
The first derivative is given by the operator, 


dx.f\ 


i,j,k 


fi—2,j,k &fi—l,j,k T 8fi+l t j,k fi+2,j,k 
12 dx 


( 2 ) 


while the second derivative is given by, 


dxxfi,j,k 


~fi+2,j,k + 16/*+1j',A; 30 fi : j t k T fi—2,j,k 

12 dx 2 ' 


(3) 
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The cross derivative is given by sequential application of 1 In the case of the xy 
derivative, this results in, 


9 X yfi,j,k 


dyfi—2,j,k 8(9 yfi—lj^k T 8dyfi + ij t k dyfi+2,j,k 

12 dx 


(4) 


We use these centered stencils for all derivatives, except advection terms for which we 
use the lop sided formulas, 


d x f, 


i,j,k 


~fi—3,j,k T 6,/i—2.j,fc l,j,k T 10 fij,k 4“ ^fi+l,j,k 

12 dx 


,6 X < 0 ( 5 ) 


&xfi,j,k 


fi+3,j,k ^fi+2,j,k T 18/i+l,j,fc 10 fi,j,k ^fi—l,j,k 

12 dx 


p x > 0 


( 6 ) 


Time integration is carried out using a fourth order accurate Runge Kutta scheme 
through the method of lines framework. A Runge Kutta step from y n to y n +\ is ac¬ 
complished with, 

Vn+l = Un + ^{kl + 2fa + 2fa + fa ), ( 7 ) 


fa = hf (x n , Y\) 
fa = hf {x n + \h, Y 2 ) 
fa = hf (x n + \h, y 3 ) 

fa = hf (x n + h, I4) 
where the quantities Y t defined by, 

Y\ = y ni Y 2 = y n + \ki, Y 3 = y n + ±k 2 , and Y 4 = y n fa fa. (9) 

serve to store boundary data. The control algorithm that is responsible for the evo¬ 
lution of the entire grid hierarchy is orchestrated by a recursive procedure, shown in 
Algorithm [I] 


Algorithm 1: A simple illustration of the AMR integration algorithm with refinement factor r. 

Procedure propagate() 

Input: int level 
Output: void 
advanceLevel(level) ; 
if level < max-level then 

foreach iteration j = 1, r do 
| propagate (level + 1) ; 

return; 
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2.4 Boundary Treatment 


The base grid consist of a ghost zone that serve to store boundary data. For refined 
grids, in addition to ghost zones, there is also a transition zone that ensures a smooth 
transition from the well resolved solution of the finer grids to the less resolved solution 
of the coarser grids. We discuss each in turn. 

2.^.1 Ghost zone 

For a fourth order accurate computation of centered derivatives, we use two ghost 
zones on each side of a gric0 Outer boundary conditions on the base grid are supplied 
by the user, for example, the user may specify periodic or sommerfeld type boundary 
conditions. We distinguish between two types of boundary conditions for the finer 
levels, those that coincide with the outer boundary and those that simply border a 
cell from the underlying coarse grid, which we refer to as refinement boundaries. If 
the ghost zones of refined grids coincide with those of the base grid, the boundary is 
filled using the prescribed procedure for outer boundaries. Otherwise, we use coarse 
grid data to fill the fine grid boundaries. 

Imposing boundary data on refinement boundaries can be reduced to the problem 
of imposing time dependent boundary conditions within a method of lines framework. 
An interesting peculiarity of higher order Runge Kutta schemes is that under certain 
conditions, they have a tendency to behave like lower order schemes, a phenomenon 
which is refered to as order reduction. This phenomenon is caused by the stiffness of 
the system under consideration in the case of implicit schemes [2B, T5 j, or inconsistent 
application of boundary data for the intermediate Runge Kutta stages in the case of 
explicit schemes [HlElj- Of course, in both cases the problem is exacerbated by the 
non linerity of the problem. Although not identified as such, order reduction effects 
in the context of mesh refinement were reported in [29IIT9] when the conventional 
method of imposing boundary conditions was used, i.e, simply applying boundaries 
corresponding to the solution evaluated at the intermediate times of the Runge Kutta 
stages. In this work, we instead use the Runge Kutta method itself to fill the ghost 
zones. We use Equations [9j where the ki are given by, [22] 

h = hy , 

h = hy' + y/ + y - f v y") , (10) 

*3 = hy' + y y" + y (j/" + S,y") . 

In the equations above, y',y" and y"' are time derivatives of the quantities under 

evolution while f y is the Jacobian matrix of the PDE system. See [22,21] for details on 

the derivation. The time derivatives can easily be obtained by polynomial interpolation 

2 However, our implementation is such that one needs the number of Ghost points to be odd for a staggered mesh and even 
otherwise. We will employ the staggered mesh in this paper and thus choose the number of Ghost points to be three. 
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methods since the coarse grid points at the advanced time will already have been 
computed before advancing the refined levels. However, a subtle issue arises in this 
case. To evolve the finer grids, at least four past points of the coarse grid solution 
are needed in order to obtain third order interpolants. This means each finer grid can 
only be initialed after the coarser grid has evolved at least four time steps. This is 
undesirable in the context of FMR. We opt to use the fact that the classical Runge 
Kutta method has a built-in interpolant, termed dense output. This interpolant is 
given by, flfi] 

4 

y(t n + Oh) =y n + E bi(0)ki + 0(/r 4 ) (11) 

2—1 

with 0 < 0 < 1 and the bi are polynomials in 0, 

bi(e)=9-^e 2 + p 3 , b 2 (9) = b 3 (e) = 0 2 - |o 3 , and b 4 (e) = -\e 2 + \e 3 

( 12 ) 

One can verify that this dense output formula reduces to Equation [7] when 0 = 1. The 
required time derivatives are then computed from Equation [Tl] as, 

§^y(tn + eh) = T £ + o(ft 4 - m ) (13) 


with 1 < m < 3. One does not need to compute the Jacobian matrix f y explicitly since 
we are only interested in the product f y y" which can be computed from the system 
TO] as, 

fyy" = ^3 ( fc 3 - fa) (14) 


See also, [21.51 S0 J. The implication is that we do not store the solution history of 
coarser grids as is currently done in some implementations. Instead, we store the four 
intermediate hi values corresponding to the current time, for all levels. This configu¬ 
ration allows for time interpolation through (jTT[) . Of course this is followed by spatial 
interpolation, for which we employ fourth order barycentric Lagrange interpolation; 
higher than fourth order was found to be unreliable in some of the runs. 


2-4-2 Transition zone 


To complete the specifications on treatment of the boundary, we examine what hap¬ 
pens close to the refinement boundary. Consider a grid hierarchy with two levels Iq 
and l\. Parametrize the solution F(x) on such a hierarchy as, 

F ( x ) = (1 -w)f(x,lo) + wf(x,h), (15) 

where f(x,lo) and f(x,h) are the solutions on the base and refined grids respectively, 
and w is a binary weight function which takes the value w = 1 if x is within the refined 
region and w = 0 everywhere else. A plot of the weight function is depicted in Figure [2j 
Note the discontinuity at the transition points x = 10 where w transitions from w = 0 





to w = 1 indicating a switch from the solution F(x) = f(x,lo) to F(x) = f(x,l i). Also 
at x = 90, w transitions from w = 1 to w = 0, indicating the switch from F(x) = f(x, l\) 
back to F(x) = f(x,lo). 


2 . 0 ,--------- 
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Fig. 2 Step function transition profile from w = 0 to w = 1. The refined region in this case is x S [10, 90]. Note the discontinuity 
at x = 10 where the weight w transitions from w=0 to w=l and again at x = 90 where w transitions from w = 1 to w = 0. 
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Because of the dispersion relation for propagating waves (17.1), there is a difference 
in phase speeds of propagating modes in the coarse and fine grid levels. As a result of 
the discontinuous transition in the weight function w, waves propagating from refined 
regions abruptly change their phase velocities when crossing refinement boundaries, 
creating a glitch that will seed spurious reflections. To circumvent this problem, we 
introduce a transition zone on the refined levels, within which the weight function w = 
w(x) is allowed to vary smoothly from w = 0 to w = 1 across the refinement boundary. 
This can be accomplished by Hcrmite interpolation. For a transition beginning at x = a 
and ending at x m b, one can derive the following profiles, 


w(a, b, x) m t 


(boxstep) 

(16) 

w(a, b , x) = 3 1 2 — 21 3 


(smoothstep) 

(17) 

w(a, b , x) = 10t 3 — 15i 4 + 6t 5 


(smootherstep) 

(18) 

where the variable t is defined as, 




f° 

¥=£ < 0 

0—CL 



‘= 1 

x—a ^ i 
b—a ^ 1 



1 x—a 

otherwise 




These profiles are shown in Figure [3} In this case the weight function varies contin¬ 
uously from w = 0 to w = 1, allowing the solution F(x) to vary smoothly across the 
transition zone. This can also be interpreted as a smooth acceleration and deceleration 
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of the associated phase speeds of propagating modes in the solution F(x). See [H] for 
a discussion on transition profiles. 

The procedure we follow to ensure a smooth transition from the fine grid solution 
to the coarse grid solution is as follows. 


- Fill ghost zone points using the method outlined in 12.4.1 


- Fill the transition zone by blending values from the refined grid, with that of the 
coarse grid according to the weight function w(x). 


The transition zone is evolved along with the fine grid solution to ensure a smooth 
coupling with the refinement boundary and thus the coarse grid solution. However, 
for reasons of stability, we do not use transition zone values when updating the coarse 
grid solution with the fine grid solution. Unless otherwise specified, we take the width 
of the transition zone to be three through out this work. For this size, the smooth 
profiles given above are equivalent. 



Fig. 3 Smooth transition profiles from w = 0 to w = 1. The refined region in this case is x E [10, 90] with the shaded regions 
representing the transition zone. We have exaggerated the width of the transition zone for ease of visualization. Compare with 
Figure [l] 


3 Evolution system 

We adopt the BSSN formulation of the Einstein held equations pJTU5]. The evolution 
equations are given in terms of the variables, 

7 ij = e ^Hj (19) 

A-ij = e Kij — - r HjK ) (20) 

subject to the constraints, 7 = Dot 7 *,- = 1, A\ = 0. Additional variables F = — 7 T 
are also introduced. The evolution equations for these variables are derived from the 
ADM equations and are given by, 
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d t (t> = —ctK + p k d k cp + ^ 0 * 0 *. ( 21 ) 

dtlij = — 2aAij + (3 k dk^ij + "ji k dj/3 k + ^j k dif3 k — ijd k /3 k , (22) 

d(K = a (AijAv + ^K 2 ^J - 7 ^'DiD ja + (3 k d k K (23) 

Qtiy = a (KAij - 2Ai k A k j) + e~ A<t, (aR ij - D,D ; n) + 

0*3*4- + 43,-0* + 40*0* - |4-3*0 fc , (24) 

^ = 2a (fi k A> k - ^Kj + 6i%-) - 2i«aj + 7 j *0\;fc + ^7 y 0* j* + 0^,i 

+ |r0V (25) 

The superscript [• • ■ ] TF denotes the trace-free part with respect to the physical metric 

7 y, and 

D/DjO = didjo - 4d (i ^)o; - f^d k a + 27 ^ 7 kl d k (j)dia. (26) 

The Ricci tensor R^ is now written as a sum of two pieces 

Rij = Rn + <•, ( 27 ) 

where Rf- is given by 

Rfj = 2l),l);o - 27jjD D/,.0 + H V4)/o - 47 ^D 0D/0, (28) 

Rij = - + 7 ro " (2r k m(i r i)ln + , ( 29 ) 


4 Numerical results 


For ease of exposition, we restrict to the case 0* = 0. Where refinement is used, we 
restrict to a refinement factor r = 2 . All runs employ three ghost points and, where 
appropriate, three transition points. 


4.1 Wave equation: Periodic boundaries 


In this section, we carry out evolutions of the wave equation with mesh refinement. 
The wave equation in flat Cartesian coordinates is given by, 


<9 tt 0 = d l di<fi. 


(30) 


ll 



We instead cast it in an alternative form, by introducing a new auxiliary variable 

n = dt<j>, 

d t (p = n (31) 

d t n = d i d i (\> (32) 

As initial data, we choose a sinusoidal profile 

<j)(x, y, z,t = 0) = sm(2n(x — t)) I7(x,y, z,t = 0) =—2n cos(2tt(x — t)), (33) 

with periodic boundary conditions on the domain x e [—0.5,0.5]. The region x e 
[—0.25,0.25] is refined by a factor of r = 2. Although the wave propagates essentially 
in one dimension, we evolve it using the full 3D grid with periodic boundary conditions 
for the outer boundaries. In Figure [4] we plot the solution errors for evolutions with 
resolution dx = l/25p for p = 1,2,3. The errors show fourth order convergence as 
desired. 
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Fig. 4 Scaled solution errors for <f> after 2 crossing times. The errors have been scaled with the resolution to highlight fourth 
order convergence. 
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4.2 Gauge Wave 

We now evolve a non linear gauge wave using the BSSN system. The Gauge wave test 
is characterized by a line element which results from a non linear gauge transformation 
of the flat Minkowski space time in Cartesian coordinates, resulting in 

ds 2 = —Hdt 2 + Hdx 2 + dy 2 + dz 2 , (34) 

where the function H is given as, 



H = H(x — t) = 1 + A sin 


2tt(x 
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for some constant A, ancl d is the wavelength. We evolve the above metric using the 
BSSN formulation with the Harmonic gauge condition, 

d t a = —a 2 K . (36) 

We choose the amplitude A m 0.1 and the wavelength d = 1. As in the last case, the 
simulation domain covers the range x e [—0.5 : 0.5] with refinement boundaries in 
the region x e [—0.25 : 0.25]. Evolving the gauge wave initial data with the BSSN 
formulation requires the addition of artificial dissipation to achieve stable evolutions. 
This is true even for unigrid runs, see for example [T] . Our dissipation operator takes 
the form 

d t Q -> d t Q + (—1 ) r/2 a J2 , (37) 

l 

for r th order accurate hnite difference stencils. 

In Figure [5] we plot the solution errors for evolutions with resolution dx = l/25p 
for p = 1,2,3,4. The errors scale according to fourth order convergence. The non 
linearity of this problem makes it ideal to demonstrate the effects of order reduction. 
Figure [6] shows a similar case, but with boundary data of the refined levels imposed 
in the conventional way. Notice the loss of convergence at high resolution. The loss of 
convergence at high resolutions was also identified in |29j . 



Fig. 5 Scaled solution errors for the g xx component of the gauge wave metric after 2 crossing times. The errors have been 
scaled with the resolution to highlight fourth order convergence. 


4.3 Wave equation: Gaussian pulse 

In the following test we show how our proposed algorithm handles artificial reflections 
that often arise when a propagating wave crosses a mesh refinement boundary. We are 
interested in waves propagating outward from the fine grid across mesh refinement 
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Fig. 6 Scaled solution errors for the g xx component of the gauge wave metric after 2 crossing times. The errors have been 
scaled with the resolution to highlight fourth order convergence. Note the loss of convergence at high resolution. 


boundaries into the coarser grid. We evolve the wave equation, with initial data given 
by a Gaussian pulse centered at the origin, 


4>(x,y, z,t = 0) = Aexp(— x 2 /a 2 ) 77(x, y, z, t — 0) = 0, (38) 


with a = 0.25 and A = 1. Although the wave propagates essentially in one dimension, 
we evolve it using the full 3D grid with periodic boundary conditions for the outer 
boundaries. The simulation domain covers x € [—4,4], The region x £ [—1,1] is further 
refined by a factor of r = 2. 


The solution is shown in Figure 7(a) The pulse starts initially at x = 0 with 


amplitude one and produces two pulses each with amplitude 0.5 traveling in opposite 
directions. In Figure [7(b) we show the result after the pulses have crossed refinement 
boundaries at x = ±1. When each pulse crosses a refinement boundary, spurious 
reflections are generated. These travel in a direction opposite that of the inducing 
pulse. When no transition zone is used, the spurious reflections reinforce at x = 0 
and can exceed the discretization error in amplitude. Employing a transition zone 
significantly reduces these artificial reflections. 


4.4 Teukolsky Wave 


In this section we evolve the Einstein held equations in three space dimensions. As 
initial data, we use the Teukolsky solution for a quadrupole l = 2, even parity m = 0 
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(a) 



Fig. 7 Solution of the wave equation at t = 2s. (a) A plot showing <f>. (b) Error in (j> computed from the analytic solution for 
two runs with the same resolution, with and without a transition zone. Note the artificial reflection at x = 0. 


waves |33j. The metric for the quadrupole modes is given by, 

ds 2 = — dt 2 + (1 + Af rr )dr 2 + (2 Bf r g)rdr dO + (2 Bf r f) r sin 6 dr d(f> 

+ (l + Cf$ + r 2 d.O 2 + [21A - 2C)h.irHu,0,Wd,p 

+ (l + Cf$ + Af®) r 2 sin 2 6 # 2 (39) 

The coefficients A, B ancf C are constructed via a generating function F(x) and are 
given by, 


A = 


'_p(2) 


3 PW 
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We take F(x) to be a superposition of ingoing (x 
waves, 


6F 

rp 5 


21 F(V 21F 

-7- 1 -c 


(40) 

(41) 

(42) 


t + r) and outgoing {x = t — r) 


F = Fft - r) + F 2 (t + r) 


and, 


p(. n ) 


d n F{x) 
dx n 


+ (-i r 
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d n F(x) 
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(43) 
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where we have chosen the particular case F\{x) = —F 2 \x) = Ae~ x . The angular 
functions f uv for even parity m = 0 modes are given by 



(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 



f m _ _,m 


(t><t> ~ Jee 


f$= 3sin 2 9-l 


We note that this represents time symmetric data and so Kij = 0 and K = = 0 


at the initial slice t = 0. We use an amplitude of A = 10 6 to complete the specification 


of initial data. 

Because of the symmetries of the problem, we impose mirror symmetry boundary 
conditions along the planes x = 0, y = 0 and z = 0. We thus evolve the Octant 
[0, 8] x [0, 8] x [0, 8] and use radiation boundary conditions at the outer boundary. We 
employ one refinement level and refine the cubic region [0,4] x [0,4] x [0,4]. The wave 
propagates radially outward crossing mesh refinement boundaries along the x = 4, 
y = 4 and z = 4 planes, eventually reaching the radiation boundary and leaving flat 
Minkowski spacetime. Although the Teukolsky wave is a routine problem for testing 
numerical relativistic codes, it is especially challenging for a mesh refinement code. 
This is mainly because the refined region is Cartesian, while the wave propagates 
spherically outward. As a result, the wavefront will not encounter the refinement 
boundaries at the same time. In Figures [8] we show the result at t = 8 for a run 
without a transition zone. Spurious ripples are generated when the wave initially hits 
the refinement boundary; the reflections continue to be generated until the wave has 
fully crossed the refinement boundary. These ripples are reflected toward the origin 
as expected. In Figure [9} we show a similar run using a transition zone. In this case, 
spurious reflections are significantly minimized. 

5 Concluding Remarks 

In this work, we have presented a fourth order convergent mesh refinement scheme 
that also significantly minimizes spurious reflections off refinement boundaries that 
are caused by differing levels of accuracy between two successive refinement levels. 
This is an important issue for the held of numerical relativity where the use of higher 
order finite differencing is becoming increasingly common [I251ITB'; :3S]. For these higher 
order methods, the truncation error can become so small that the dominant error 
comes from spurious reflections. Our method is not restricted to any formulation of 


16 



le-06 


-le-06 


-1 



8e-07 

6e-07 

4e-07 

2e-07 

0 

-2e-07 

-4e-07 

-6e-07 

-8e-07 


Fig. 8 Evolution of the 7 zz component of the metric along the xy plane without a Transition zone. Note the spurious ripples 
in the refinement region. For ease of visualization, we mark the boundary of the refined grid with with white lines. 
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Fig. 9 Evolution of the 7 component of the metric along the xy plane with a Transition zone. Note the absence of spurious 
ripples in the refinement region. For ease of visualization, we mark the boundary of the refined grid with with white lines, 
compare with Figure |ff[ 


the Einstein field equations. Indeed one can apply it to any hyperbolic system of 
partial differential equations. 

Because we are not using the buffer zone approach as detailed in pMITj. our method 
requires a total of six points on each boundary, irrespective of the time integration 
method used. This differs from employing buffer zones in that, for a fourth order ac¬ 
curate Runge Kutta algorithm, along with fourth order finite differencing involving 


17 











lop-sided advection stencils, twelve points are needed along each boundary (the situ¬ 
ation could be worse for higher order finite differencing) 11~H. 7i. This is a significant 
saving in memory usage, especially in three space dimensions where the buffer zone 
can be a significant part of the grid. In addition, the blending operation we employ 
to fill the transition zone is cheaper than having to repopulate the entire buffer zone 
after every time step. We also note that the use of a transition zone is computation¬ 
ally cheaper than the sponge boundary method since one has to populate the sponge 
boundary at every intermediate Runge Kutta step, while the transition zone is only 
populated at the end of the time step. Moreover, there is often a level of experimenta¬ 
tion required to determine how large a sponge zone one should use. Although one can 
have the transition zone as large as desired, we have found satisfactory results with a 
size that spans only three fine grid points. 

The implementation described here uses Runge Kutta dense output formulas to 
interpolate in time, which avoids any potential issues with polynomial interpolation. 
In particular, because one does not have to couple fine grid solutions to the solution 
history of coarser grids, fine grids can be immediately initialized along with the base 
grid. This is especially attractive since fine grids can be initialized to the same ac¬ 
curacy as the base grid, by using the same initialization routine as the base grid. In 
addition, the use of high order polynomial interpolation in time requires the use of 
asymmetric stencils, which results in oscillations in the refinement boundaries and ul¬ 
timately, unstable modes at higher resolutions [7j. This dense output method of time 
interpolation, which is routinely used in the solution of ordinary differential equations 
|IB], was also used in the context of mesh refinement pTJ for conservation laws and 
j9] for solving Maxwell’s equations. 

It would be interesting to investigate the efficiency of the transition zone imple¬ 
mentation in an adaptive context or whether it would minimize spurious reflections 
that are a result of shock waves crossing refinement boundaries. There is also the 
question of how it would affect conservation along interface boundaries in the context 
of conservation laws. These and other issues involving black hole spacetimes are a 
subject of further study. 
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7 Appendix 

7.1 Dispersion relation 

In this section we calculate the dispersion relation and phase velocity resulting in 
discretizing the wave equation with a fourth order stencil along with Runge Kutta time 
marching. We follow the approach (and notation) of jTO] where a similar calculation 
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was given using second order finite differences and iterative Crank-Nicholson time 
marching. 


The wave equation (31) can be written in matrix form as 

v ‘ = (lo) V ’ 

where we have defined the vector V as 


1/ = 


n 


(53) 


(54) 


With this identification, we denote by V 7 ] the solution at time step n and grid point 


j. The second derivative operator d xx appearing in (J53J) is given by the stencil, 

n 

2—1 


a T , n _ ~V n i + 2 + 16^ + i - 30W* + 16V n ,_, - Y \; 

OXT, V r! - 


2-2 


12 dx 2 

For our analysis, we consider plane wave solutions for V, of the form, 

yn _ yy e iundt e ~ikjdx 


(55) 


(56) 


for some constant vector W. Using the classical fourth order Runge Kutta scheme to 

(57) 


advance (53) in time, results in the update rule, 


V n+1 . = MV n 
j r 


Plugging in (56) to the above rule results in the relation, 

1 — 2T 2 + |d 4 dt(l~lA 2 ) 

3dt 


( 1 - 2d 2 + §A 4 dt (1 - 
’" d ‘ w = ( 4/l 2 (3 - 2/1 2 ) j _ 2a2 

\ 3 dt 


+p 4 


w 


(58) 


where we have defined A as, 


s\W in 2 (^)-n sin2(w *) 


(59) 


The system (58) represents an eigenvalue problem. In particular, W is an eigenvector 
corresponding to the eigenvalue e luldt for the matrix in (58). Further analysis shows 
the eigenvalues to be the pair, 


: ,ljJi = 1-2 A 2 + |/1 4 ± 2iA ( 1 - |yl 2 


(60) 


This expression represents the dispersion relation, relating the frequency c o with the 
wave number k. For completeness we calculate the phase velocity, v p (X) = i/k for 


£ = Re (cc). From (60) we get, 


£dt = arcsin 


2/1(3 - 2 A 2 ) 
^9 - 4T 6 (2 - A 2 ) 


(61) 
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Therefore, 


v p (X) 


X 

27 rdt 


arcsin 


2/1(3 - 2 A 2 ) 
^9 - 4/«(2 - A 2 ) 


(62) 
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